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Abstract: This paper proposes a radial dependent dispersive finite- 
difference time-domain method for the modelling of electromagnetic 
cloaking structures. The permittivity and permeability of the cloak are 
mapped to the Drude dispersion model and taken into account in dispersive 
FDTD simulations. Numerical simulations demonstrate that under ideal 
conditions, objects placed inside the cloak are 'invisible' to external 
electromagnetic fields. However for the simplified cloak based on linear 
transformations, the back scattering has a similar level to the case of a PEC 
cylinder without any cloak, rendering the object still being 'visible'. It is 
also demonstrated numerically that the simplified cloak based on high-order 
transformations can indeed improve the cloaking performance. 
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1. Introduction 

The widespread interest in the invisibility of objects has led to the recent development in elec- 
tromagnetic cloaking structures. Pendry et ah proposed an electromagnetic material through 
which electromagnetic fields can be controlled and manipulated to propagate around its inte- 
rior region like the flow of water [1], hence objects placed inside become 'invisible' to ex- 
ternal electromagnetic fields. The proposed cloaking structure in (TJ requires inhomogeneous 
and anisotropic media, with both permittivity and permeability independently controlled and 
radially dependent. The magnitudes of the relative permittivity and permeability of the perfect 
cloak are less than one, therefore it cannot be constructed using naturally existing materials. The 



recent development of metamaterials 1 2 ] (artificially engineered structures with exotic electro- 
magnetic properties which cannot be obtained naturally) allows the control of material parame- 
ters with freedom, and hence the construction of such cloaking structures. However same as the 
negative-index metamaterials [2], the cloaking materials are inevitably dispersive and therefore 
band-limited. Furthermore, the complete set of material parameters proposed in [1 j requires the 
control of all the components of permittivity and permeability of the material, making it dif- 
ficult for practical realisation. It has been proposed to use reduced sets of material parameters 
for both transverse electric (TE) [3] and transverse magnetic (TM) (4) cases. The reduced set 
of material parameters for the TM case eliminates the dependence on the material's magnetic 
property, which is especially important for the realisation of the cloak in the optical frequency 
range due to the absence of optical magnetism in nature. However, considerable reflections oc- 
cur because of the impedance mismatch at the outer boundary of such a simplified cloak. Under 
the assumption of the geometric optics, a high-order transformation has been proposed in 
to improve the performance and minimise the scattering introduced by the cloak. 

The development of Pendry's cloak is based on the coordinate transformation technique (U 
[6) , which has also evoked other research topics such as the design of magnifying perfect and 
super lenses J7), the transformation media that rotate electromagnetic fields [8], the design 
of reflectionless complex media for shifting and splitting optical beams |9], and the design 
of conformal antennas ifTOl etc. The spatial transformation technique has also been applied to 
anslyse eccentric elliptic cloaks ITU and for acoustic cloaking |[T2lfT3ll . Other theoretical studies 
of the cloaking structure include the sensitivity of an ideal cloak to tiny perturbations llT4l . the 
performance of the cylindrical cloaks with simplified material parameters [1131 [T6l [TTl . the 
realisation of cloaking using a concentric layered structure of homogeneous isotropic materials 
|[T8l . the improvement of the cloaking performance using soft-and-hard surface lining lfT9l , and 
the broadband cloaking using sensors and active sources near the surface of a region l20l etc. 
The experimental demonstration of a simplified cloak consisting of split-ring resonators (SRRs) 
has been reported in microwave frequencies I2H . For the optical frequency range, the cloak can 
be constructed by embedding silver wires in a dielectric medium (4), or using a gold-dielectric 
concentric layered structure li22l . 

It is worth mentioning that there exist different approaches to render objects invisible, for 
example, by canceling the dipolar scattering using plasmonic coatings (231 [24l, and using a 
left-handed media (LHM) coating ll25ll26l . However, the plasmonic coating approach is limited 
to the sub- wavelength scale of the object, and the coating depends on the geometry and material 
parameters of the object to be cloaked; the cloaking performance using LHM coating is affected 
by the objects placed inside. In comparison to these different approaches, Pendry's cloak is 
more general and can be applied to objects with any dimensions and under any wavelength 
condition. 

The modelling of Pendry's invisible cloak has been performed by using both analyti- 
cal and numerical methods. Besides the widely used coordinate transformation technique 
d [6j HI [U HI HOl EH H2l H3l HSl Ell, a cylindrical wave expansion technique d, and a 
method based on the full- wave Mie scattering model ll28l have also been applied. In addition, 
the full- wave finite element method (FEM) based commercial simulation package COMSOL 
Multiphysics™ has been extensively used to model different cloaks and validate theoretical 
predictions H B El ESI due to its ability of dealing with anisotropic and radial depen- 
dent material parameters. While most of the numerical simulations have been performed in 
the frequency domain, little attention has been paid to the time-domain analysis of the cloaks. 
Frequency domain techniques such as FEM may become inefficient when wideband solutions 
are needed. So far the only time-domain analysis of the cloak is presented in ll29l using a 
time-dependent scattering theory. In this paper, we propose a dispersive finite-difference time- 



domain (FDTD) method to deal with both the frequency and radial dependent permittivity 
and permeability of the cloak. Due to its simplicity in implementation and ability of treating 
anisotropic and inhomogeneous materials, the FDTD method has been extremely popular for 
the analysis of waveguides, transmission lines and antennas as well as the calculation of disper- 
sion diagrams of photonic crystals in both microwave and optical communities. The aim of this 
paper is to provide a full- wave FDTD solution to assist the analysis, design and optimisation of 
the cloaking structures. To the authors' knowledge, it is the first time that the FDTD method is 
implemented to model cloaking structures. 

2. Dispersive FDTD modelling of the cloaking structure 

A complete set of material parameters of the ideal cloak is given by Q : 

r-Ri r ( R 2 \ 2 r-R x 

e r = Hr=— 8^=^ = —, * = ft =(^J — -, (1) 

where R\ and R 2 are the inner and outer radius of the cloak, respectively. It can be easily 
identified from (Q]) that the ranges of permittivity and the permeability within the cloak are 
£ r ^ r e[0,(R 2 -R l )/R 2 ],£^^ e [R 2 /(R2-Ri)M and £ z ,\i z e [0,R 2 / (R 2 -Ri)]- Since the 
values of £ r , /i r , £ z and \i z are less than one, same as the case for the LHMs, the cloak cannot 
be directly modelled using the conventional FDTD method with constant material parameters 
(at one particular location r). However, one can map the material parameters using dispersive 
material models e.g. Drude model (taking £ r as an example) 

col 

£r((0) = l- P - , (2) 

co 2 -jcoy 

where co p and y are the plasma and collision frequencies of the material, respectively. By vary- 
ing the plasma frequency, the radial dependent material parameters (Q} can be achieved. Note 
that in practise, the plasma frequency of the material depends on the periodicity of the SRRs 
l2T1l or wires EI, which varies along the radial direction. Furthermore, different dispersion 
models (e.g. Debye, Lorentz etc.) can be also considered for the modelling of the cloak, which 
will lead to slightly different FDTD formula from the following ones. 

Since the conventional FDTD method (30l [31]] deals with frequency-independent materials, 
the frequency-dependent FDTD method is hence referred as the dispersive FDTD method l32l 
[33l [341 . For simplicity, in this paper, we have implemented the dispersive FDTD method for 
a two-dimensional (2-D) case and only three field components are non-zero: E x , E y and H z . 
The modelled cloak is cylindrical and infinite along z-direction (fi r = fi^ = £ z = in (Q])), 
however an extension to a three-dimensional (3-D) FDTD method to model a 3-D cloak (TJ is 
straightforward. There are also different dispersive FDTD methods using different approaches 
to deal with frequency-dependent material parameters: the recursive convolution (RC) method 
(32), the auxiliary differential equation (ADE) method ll33l and the Z- transform method l34l . 
Due to its simplicity, we have chosen the ADE method for the modelling of the cloak. 

The ADE dispersive FDTD method is based on Faraday's and Ampere's Laws: 

VxE = (3) 
VxH = ^ (4) 



as well as the constitutive relations D = £E and B = /iH where £ and \i are expressed by ©• 
Equations © and © can be discretised following a normal procedure (30l [5U which leads to 



the conventional FDTD updating equations: 

B n+1 



= B* 
= D n 



■A/-VxE" + 2, 



(5) 
(6) 



where V is the discrete curl operator, At is the discrete FDTD time step and n is the number of 
the time steps. 

In addition, auxiliary differential equations have to be taken into account and they can be 
discretised through the following steps. For the conventional Cartesian FDTD mesh, since the 
material parameters given in (Q]) are in cylindrical coordinates, the coordinate transformation 



£ r cos 2 + £</> sin 2 (j) [e r — £</>) sin (j) cos (j) 



is used. The tensor form of the constitutive relation is given by 
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Note that the inverse of the permittivity tensor matrix exists only when £ r ^ and £^ ^ 0, 
which is not the case for the inner boundary of the cloak. In our FDTD simulations, we place a 
perfect electric conductor (PEC) cylinder with radius equal to R\ inside the cloak to guarantee 
the validity of ©. 

Substituting © into (8]) gives 



(10) 



£ r £0 £ E X = (e r sin 2 (j) + £0 cos 2 (j)) D x + - £ r ) sin cos 0/^ 

£ r £0 £b£y = (fir COS 2 (j) + £^ sin 2 0) + (£0 — £ r ) SUl COS 0Z) X ' 

Express £ r in the Drude form of (Q]), Eq. dTOb can be written as 

£ O £0 (co 2 - jcoy- co 2 ) £ x = [(co 2 - jcoy- co 2 ) sin 2 + £</> (^ 2 - yffly) cos 2 0] £> x 

+ [£0 (co 2 - jcoy) - (co 2 - jcoy- co 2 )] sin cos (j)D y , 
£ O £0 (a) 2 - jcoy- co 2 ) E y = [(co 2 - jcoy- co 2 ) cos 2 + £$ (co 2 - jcoy) sin 2 0] 

+ [£ (co 2 -jcoy) - (co 2 - jcoy- co 2 )] sin0cos0£> x . 

Notice that £^> remains in (fTTb because its value is always greater than one (except at the inner 
surface of the cloak) and can be directly used in conventional FDTD updating equations (30l 
|3H . Using an inverse Fourier transform and the following rules: 



(11) 



j co 



d 2 d 2 

the first equation of (TTTh can be rewritten in the time domain as 



(9 2 d 
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sin (j) cos (j>D y . (13) 



The FDTD simulation domain is represented by an equally spaced 3-D grid with periods 
Ax, Ay and Az along x-, y- and z-directions, respectively. For discretisation of Eq. (TDK we use 
the central finite difference operators in time (8 t and 8 2 ) and the central average operator with 
respect to time {ji t and nf): 
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where the operators 5t, 8 2 , jU, and fi? are defined as in ll35l : 
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Here F represents field components and m Xj m y ,m z are indices corresponding to a certain dis- 
cretisation point in the FDTD domain. The discretised Eq. (IT3l) reads 
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Note that in ([T5lh the discretisation of the term co 2 of (TT31) is performed using the central average 
operator /I 2 in order to guarantee an improved stability; the central average operator \i t is used 
for the term containing y to preserve the second-order feature of the equation. Equation (TT3T) 
can be written as 
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Therefore the updating equation for E x can be obtained as 
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where the coefficients a x to l x are given by 

n2l 



sin 2 (j) 



sin 2 (j) 



sin 2 (j) 



(At) 2 2At 4 



2 . <°p 
(AO 2 2 



_J r_, % 

(At) 2 2At 4 



+ 8(f) COS 2 


1 

JKtf 


2Ar_ 












+ £(j, COS 2 


l 


7 " 
2Ar_ 



£0 
£0 



^_+x 

(At) 2 2At 



2 " 




2 


(A?) 2 . 




(AO 2 



(AO 2 2Ar 4 



_J 7_ 

(Ar) 2 2Af 



_J L , _i 

(At) 2 2At 4 



"(Ar) 2+ 2 



, h x = e e^ 



1 



sin cos , 
), 

sin cos , 
7 



(At) 2 2At 4 

Following the same procedure, the updating equation for E y is 
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where the coefficients <z-y to are given by 
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Note that the field quantities D x in (TT8t and Z>y in (fTTI) are locally averaged values of D x and Z>y, 
respectively since the x- and y-components of the field are in different locations in the FDTD 
domain. The averaged values can be calculated using 



D x (i,j) 



D x (i, j) + Djj, j + 1 ) + D x {i -l,j)+D x (i-l,j+l) 
4 

D y (iJ)+D y (i+ \,j)+D y (i,j-\)+D y {i+\J-\) 



(19) 



where (i, j) is the coordinate of the field component. 



The magnetic field H z can be calculated starting from the constitutive relation B z = fl z H z , 
where jl z is defined in (Q]). Again we map the permeability using the Drude model 



ji z {co)=A \ 1 



CO 



pm 



(20) 



where we choose A = 27^/(^2 — R\), and co pm and y m are the magnetic plasma and collision 
frequencies of the material, respectively. Following the above discretisation procedure, the up- 
dating equation for H z can be derived as 




.(21) 



Equations ©, ©, (fTTlh (IT8l) and (I2T1) form an FDTD updating equation set for the well-known 
leap-frog scheme l30l . 

Since the FDTD method is inherently a numerical technique, the spatial as well as time 
discretisations have important effects on the accuracy of simulation results. Also because the 
permittivity /permeability is frequency dependent, one can expect a slight difference between 
the analytical and numerical material parameters due to the discrete time step in FDTD. In gen- 
eral, for the modelling of conventional dielectrics with the relative permittivity/permeability 
greater than one, a spatial resolution (FDTD cell size) of Ax < A/ 10 is required [31 1. However 
from our previous analysis |[38l that for the modelling of metamaterials especially LHMs, the 
numerical errors due to the time discretisation will cause spurious resonances, hence a require- 
ment of Ax < A /80 is necessary. For the case of the cloak, following the same approach as in 
and for the case of plane- waves, substituting 



E n = Ee jncoAt , 



-T>e jn(oAt , 



(22) 



into Eq. (fT6K and comparing the resulting equation with the first equation of (ITOb one can 
find that (numerical permittivity along (j) direction) has the exact analytical value, while e r 
(numerical permittivity along radial directions) takes the following form: 

^ _ ( co 2 p (At) 2 cos 2 ^ 1 

8r ~ e ° 1 ~ 2 sin ^ (2 sin ^ - jyAtcos ~^f) ' ^ 

Note that Eq. (f23l) simplifies to the Drude dispersion model Q when At — > 0. We plot the 
comparison between the analytical © and numerical relative permittivity (f23l) for the case of 
£ r = 0.1 (lossless) in Fig.Q] It is apparent that the conventional requirement of FDTD spatial 
resolution does not guarantee the accuracy for the modelling of the cloaks, and even for the 
case of Ax = X /40, the numerical error is still about 2%. 

With the expression of the numerical permittivity (l23l) available, one can correct the errors 
introduced by the discrepancy between numerical and analytical material parameters. For ex- 
ample, if the required permittivity is £ r = e' r + j£ f r f , after simple derivations, the corrected plasma 
and collision frequencies can be calculated as 
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Fig. 1. The comparison between the analytical (2]) and numerical material parameters (fzTb 
for different FDTD spatial resolutions for the case of e r = 0.1 (lossless). 



It is found from our FDTD simulations with Ax = A/35 that without the correction of nu- 
merical material parameters, the simulation becomes unstable before reaching the steady- state. 
The cause of such an instability is currently under our investigation. On the other hand, the 
correction of material parameters (both permittivity and permeability) guarantees stable FDTD 
simulations. Therefore in the following, the corrected material parameters (l24t are always used. 

3. Numerical results and discussion 

The 2-D FDTD simulation domain is shown in Fig. [2 The following parameters are used in 
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Fig. 2. A two-dimensional (2-D) FDTD simulation domain for the case of plane-wave in- 
cidence on the cloak. 

the simulations: FDTD cell size Ax = Ay = A/ 150 where A is the wavelength at the operating 
frequency / = 2.0 GHz, and the time step At = Ax/ V2c, chosen according to the Courant sta- 
bility criterion [31]. In the present paper, we assume the ideal lossless case i.e. the collision 
frequency in © is equal to zero (7=0). The radial dependent plasma frequency can be calcu- 
lated using (l24l) with a given value of £ r calculated from (Q]). The radii of the cloak are: R \ = 0. 1 
m and R2 = 0.2 m. The computation domain is truncated using Berenger's perfectly matched 
layer (PML) l37l in y-direction to absorb waves leaving the computation domain without intro- 



ducing reflections, and terminated with the periodic boundary conditions (PBCs) in x-direction 
for the modelling of a plane- wave source. The plane- wave source is implemented by specifying 
a complete column of FDTD cells using a certain wave function (sinusoidal source in our case), 
see Fig. [2 

First we consider the ideal cloak, whose material parameters are given by with \i r = 
M</> = e z — 0- Figure [3] shows the distributions of the electric and magnetic field components 
calculated from the dispersive FDTD simulation of the ideal cloak. Note that only the central 
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Fig. 3. (a) Material parameters for an infinite ideal cylindrical cloak [3] where all £ r , 
and ji z are radial dependent, (b), (c) and (d) Field distributions from dispersive FDTD sim- 
ulations of the cloak: (b) x-component of the electric field, (c) v-component of the electric 
field and (d) the magnetic field. 

part of the simulation domain is shown and the actual computation domain is larger. It is seen 
in Fig. Ob), (c) and (d) that the distribution is in good agreement with those in (3l|4l calculated 
using frequency domain methods, although there are tiny ripples in the magnetic field. In fact 
the field distribution presented in shows even stronger ripples which may be due to the 
inadequate spatial resolution used in the calculation. These ripples are purely numerical and 
will disappear when an extremely fine FDTD grid is used. In our case, the slight ripples are 
also caused by the staircase approximation of the circular surface of the cloak in the Cartesian 
coordinate system. There are also non-zero scattered fields in the x-component of the electric 
field outside the cloak (see Fig. Ob)), which is also a sign of numerical errors because an ideal 
cloak does not introduce any scattering outside the cloak. The staircase approximation only 
causes a very small amount of numerical errors due to the fine mesh used in our simulations. 
The staircase problem can be further alleviated to improve the accuracy of simulations using a 
conformal scheme in addition to the dispersive FDTD method i.e. using an effective permittivity 
at material boundaries, as it was done for the case of isotropic dispersive materials at planar 



and curved boundaries [ 39 1 . However, the difficulty lies in the anisotropy of the cloaking 
material, which leads to an eighth-order differential equation to be discretised (the order of the 
differential equation for the case of isotropic dispersive materials is four ll39l ). 

As mentioned previously, it is difficult to realise the ideal cloak with all the components of 
permittivity and permeability being radial dependent. Therefore it is proposed in [4] that while 
keeping the same wave trajectory, a reduced set of material parameters also allows to construct 
a simplified cloak. The reduced set of material parameters is given by H : 



R2—R1 



£(j) 



R2 —R\ 



1. 



(25) 



Following the same procedure, we have also modified the above proposed dispersive FDTD 
method to model the reduced set of material parameters and analysed its cloaking performance. 
The dimensions of the simplified cloak remain the same as the ideal one. The field distributions 
in the steady- state of the simulation are plotted in Fig. HI Such a cloak significantly reduces 
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Fig. 4. (a) Material parameters for an infinite simplified cylindrical cloak using a linear 
transformation J4) where only e r is radial dependent, (b), (c) and (d) Field distributions 
from dispersive FDTD simulations of the cloak: (b) x-component of the electric field, (c) 
v-component of the electric field and (d) the magnetic field. 



the complexity of practical realisation since only e r is radial dependent, as shown in Fig. |4^a). 
However, considerably reflections occur due to the fact that the impedance matching at the outer 
boundary of the simplified cloak is not satisfied anymore, see Fig. HJb) and (d). Interestingly, 
the v-component of the electric field is only affected slightly by the scattered field. Note that 
here we only consider the simplified non-magnetic cloak Q], the simplified cloak proposed in 
||3l can be modelled in a similar way. 



The scattering from the simplified cloak (l25t due to the impedance mismatch can be reduced 
by using an improved cloak based on a high-order transformation [5]. The material parameters 
are given by 



£0 



dg{r>) 



1 -2 



(26) 



where r = g(r') — [(R\/R2)(r J /R2 — 2) + 1] r' -\-R\. We have also modelled such a cloak and 
plotted the field distributions in Fig. [5] Its dimensions are kept the same as the previous two 
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Fig. 5. (a) Material parameters for an infinite simplified cylindrical cloak using a high-order 
transformation Q where only e r and are radial dependent, (b), (c) and (d) Field distribu- 
tions from dispersive FDTD simulations of the cloak: (b) x-component of the electric field, 
(c) y-component of the electric field and (d) the magnetic field. 

cases. In fact the dimensions of this cloak is at its limit since it is required to have R\ /R2 < 0.5 to 
guarantee a monotonic transformation [5) . The improved cloak adds an additional dependency 
of the permittivity on the radius, as shown in Fig. Oa). It is clear that indeed the high-order 
transformation based cloak has an improved impedance on its outer boundary and hence con- 
siderably reduces the scattered field. Notice that the wavefront only starts to bend near the inner 
surface of the cloak, in comparison to the case for the ideal cloak shown in Fig. [3] This is due 
to the slow variance of the impedance from the outer boundary towards inside the cloak. 

For demonstration, we have also plotted the power flow diagrams as well as the scattering 
patterns in Fig.[6]for the above ideal cloak, the simplified cloak based on a linear transformation, 
and the simplified cloak based on a the high-order transformation. The power flow diagrams 
show that for the ideal case (Fig. 0a)), the wavefront enters the cloak smoothly, bends around 
the central region and returns to its original pattern after leaving the cloak. The cloak based 
on the high-order transformation (Fig. 0c)) shows a similar pattern with a smooth bending of 
the wavefront near the central region of the cloak. However due to reflections, the power flow 
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Fig. 6. Power flow diagrams for (a) the ideal cloak Q, (b) the simplified cloak based on the 
linear transformation J4) and (c) the simplified cloak based on a high-order transformation 
l5l. (d) Comparison of the scattering patterns for different cloaks and for the case of the 
PEC cylinder without cloak: black solid line - PEC cylinder; yellow dot-dashed line - linear 
transformation based simplified cloak; blue dashed line - high-order transformation based 
simplified cloak; red solid line - ideal cloak. 

is disturbed before entering the cloak, while the wavefront leaving the cloak has a relatively 
smoother distribution. For the case of the cloak based on the linear transformation as shown in 
Fig. Ob), the waves do not strictly follow the trend of the bending inside the cloak and propagate 
in arbitrary directions. The consequence is that the external field is significantly disturbed. This 
can be clearly identified from the scattering patterns. The scattering pattern is calculated with 
the reference to the free space case without any cloak or objects, and then normalised to the 
scattering pattern of a PEC cylinder without any cloak, see Fig. [Ud). For all the cloaks, the 
scattering at the back of the cloak (relative to the direction of wave incidence) is dramatically 
reduced. However the level of scattered field for the linear cloak is almost the same as the case 
for a PEC cylinder without any cloak, leading to the conclusion that the object placed inside 
this simplified cloak can be detected from its front side, similar to the conclusion drawn in 
[1131 . For the high-order transformation based cloak, the scattering is reduced by around four 
times comparing with the linear transformation based one. Theoretically the ideal cloak has 
zero scattered field, however the nonzero (small) values in Fig. Od) is mainly caused by the 
staircase approximation in FDTD simulations, and will approach zero when an extremely fine 
FDTD grid is used or a conformal scheme is employed, as mentioned earlier. 



4. Conclusion 



In conclusion, we have proposed a dispersive FDTD method for the modelling of the cloaking 
structure. The unusual material parameters (the relative magnitudes of the permittivity and per- 
meability are less than one) are mapped to the Drude dispersion model, which is then taken into 
account in FDTD simulations using an auxiliary differential equation based method. The pro- 
posed method is implemented in a two-dimensional case and three different cylindrical cloaks 
are considered in our simulations: the ideal cloak, the linear transformation based cloak and the 
high-order transformation based cloak. It is found from the simulations that the linear transfor- 
mation based cloak introduces a level of back scattering similar to the one of a PEC cylinder 
without the cloak, causing the possibility of the object being detected. Such scattering can 
be significantly reduced by using the high-order transformation based cloak. In this paper we 
have only considered lossless cloaks. The 'ideal' cloak with a material loss of tan 5 = 0.1 has 
been modelled in using the finite element method, and the case for tan 5 =0.01, 0.1 and 
1 modelled in [28 ] using a full- wave Mie scattering theory. Lossy cloaks can be also directly 
modelled using the above proposed dispersive FDTD method by specifying a certain value for 
the collision frequency in the Drude model for £ r and defining a dielectric loss for . 



